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Open quantum systems that interact with structured reservoirs exhibit non-Markovian dynamics. 
We present a quantum jump method for treating the dynamics of such systems. This approach is a 
generalization of the standard Monte Carlo Wave Function (MCWF) method for Markovian dynam- 
ics. The MCWF method identifies decay rates with jump probabilities and fails for non-Markovian 
systems where the time-dependent rates become temporarily negative. Our non-Markovian quan- 
tum jump (NMQJ) approach circumvents this problem and provides an efficient unravelling of the 
ensemble dynamics. 
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Introduction. When an open quantum system inter- 
acts with a reservoir having non-trivial structure, the 
system dynamics exhibits non-Markovian memory ef- 
fects [T]. The information on the state of the open system 
is contained in the density matrix whose time evolution 
is governed by a master equation consisting of two parts. 
The system Hamiltonian induces unitary evolution of the 
density matrix while the dissipative part, which includes 
the information on the properties of the environment in 
form of decay rates, induces non-unitary effects via the 
jump operators. Already for Markovian systems, which 
do not have memory, finding the solution of the mas- 
ter equation may be very complicated. The task gets 
even more challenging with non-Markovian systems and 
structured reservoirs. Such systems display modified de- 
cay dynamics paving the way to new types of quantum 
control techniques [2]. 

Non-Markovian systems appear in many branches of 
physics, such as quantum optics (HOE], solid state 
physics jlj, and quantum chemistry [5]- In quantum in- 
formation processing [6] , the non-Markovian character of 
decoherence has to be accounted for and it leads to the 
concept of non-Markovian quantum channels 7 . De- 
coherence also plays a central role in the transition from 
quantum to classical world [5]. In fact, non-Mar kovianity 
has been recently proposed as a means to manipulate 
the quantum-classical border [9J. Since it is elusive to 
solve the open system dynamics, new methods for non- 
Markovian systems are highly desired. 

In this Letter we provide an efficient way to unravel 
a general non-Markovian master equation. The different 
ways to build an ensemble of stochastic wave functions 
describing the density matrix fall roughly into two cat- 
egories [10]: time-evolution including (i) discontinuous 
changes (quantum jumps), e.g., the Monte Carlo Wave 
Function (MCWF) method [TT]; (ii) continuous stochas- 
tic changes, e.g., the Quantum State Diffusion (QSD) 
method |12L 113] - Our non-Markovian quantum jump 
(NMQJ) method generalizes the widely used Markovian 
MCWF into the field of non-Markovian systems, and thus 
belongs to the first of the two categories. 



There exists a non-Markovian variant of QSD [J2] 
and a somewhat related formulation [14 . These meth- 
ods, however, are difficult to implement beyond very 
simple examples. Other unravelings of non-Markovian 
master equations contain fictitious harmonic oscillator 
modes [TS] and pseudomodes [15], or some other forms 
of extensions of the system Hilbert space [TTJ, [TH] . One 
formulation, using quantum jumps, exploits an analogue 
to the the hidden variable theory [19]. The use of ex- 
tended Hilbert spaces comes always with an added cost 
for computational efficiency. 

Our formulation avoids the use of Hilbert space ex- 
tensions and is based on the following observation. The 
information, which the system loses to the environment 
at the time of the jump, can be later recovered by the 
system due to non-Markovian memory. We show explic- 
itly how this happens on the level of single realizations. 
Before discussing on the insight and benefits that our 
NMQJ method provides, we first introduce the master 
equation and the method, and present a case study with 
an atom in a photonic band gap. 

Non-Markovian master equation. The non-Markovian 
dynamics of the reduced system density matrix p(t) is 
given by the master equation [1] 

p(t) = ^[H s ,p{t)]+Y / ^At)c j {t)p{t)c){t) 
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- ^A,(t){p(<),C](t)Q(i)}. (1) 
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Above, H s is the system Hamiltonian and Cj(t) are the 
jump operators describing changes in the system due to 
interaction with the reservoir. Aj(t) is the decay rate of 
channel j. It can be shown that the most general master 
equations local in time for non-Markovian systems can 
be cast in the form of Eq. (JTJ 18J. In the Markovian 
case all Aj are positive constants. In the non-Markovian 
case the rates may oscillate and take negative values for 
finite time intervals. This is a sign of the non-Markovian 
memory effects and reflects the exchange of information 
back and forth between the system and the reservoir. 
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MCWF and NMQJ methods. The system proper- 
ties are calculated as an average over the state vector 
ensemble of size N and we follow closely the MCWF 
method A generic way to write the density matrix 
is 



pit) = E 



N a (t) 
N 



\Mt))(Mt)l 
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where N a (t) is the number of ensemble members in the 
state \ip a (t)) at time t. The deterministic evolution of a 
given state vector \ijj a (t)), for small enough time steps St 
and before the renormalization, is given by 



\4> a (t + 6t))= 1- 



iHSt 



(3) 



where the non-Hermitian Monte Carlo Hamiltonian is 
H = H s - ihJ2 3 l&jtyCjityC^t) and the renormal- 
ized state is \ip a (t+St)) = \<j> a {t + St))/\\\<p a (t + St))\\. 
For positive decay channels j+, Aj + (i) > 0, the deter- 
ministic evolution is interrupted by jumps \tp a (t)) — > 
Cj + (t)\t/; a (t)) /\\Cj + (t)\ip a (t))\\ which occur with proba- 
bility 

Pi+(t) = A j+ (t)5t(Mt)\Cl(t)C j+ (t)\<p a (t)), (4) 

during time step St [11] . The Markovian MCWF method 
can be extended to the situations where the rates become 
time dependent, but this is limited to positive decay rates 
only. 

In our approach the non-Markovian quantum jumps 
for negative channels Aj_ (t) < 0, have the form 
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(5) 



where the source state of the jump is \ip a {t)) — 
Cj_ {t)\i) a ,(t)) /\\C _{t)\i) a ,(t))\\. This transition for a 
given state vector \ip a ) in the ensemble ([2| occurs with 
the probability 
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Note that the probability of the non-Markovian jump is 
given by the target state of the jump. The sign of the 
decay rate Aj (t) can be understood in the following way. 
First, when for a given channel j, Aj(t) > 0, the process 



goes as \ip c 



\1>o 



QhM/IIQI^o 



Later on, 



when the decay rate becomes negative, Aj(t) < 0, the 
direction of this process is reversed and the jump occurs 
to opposite direction |f/v) iV'a)- 

The proof of our NMQJ method goes in a very simi- 
lar way to that of the Markovian MCWF method [IT] . 
By weighting the deterministic and jump paths over time 



step St with the appropriate probabilities we should ob- 
tain the master equation ([lj. Calculating the average a 
of the evolution of the ensemble (l2| over St gives 
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Here, the summations a and a' run over the ensemble, 
the summation over j + and j'_ cover the positive and neg- 
ative channels, respectively. The first term on the r.h.s. 
for the summation over a is the product of the no-jump 
probability and the deterministic evolution of the state 
vector, the second and third terms describe the positive 
and negative channel jumps, respectively, with the cor- 
responding probabilities. By using Eq. the last term 
can also be written as J2j_, a > ^a^a'WlVvWK^a'WI- 
In general, using Eqs. in Eq. ^fj, and keeping in 

mind the form of the reversed jump \ij) a ) — ► \ip a ') with 
l^a) — Cj_ \ip a ')/\ \4>a') II, gives the master equation 

Example: Photonic band gap. To illustrate the NMQJ 
method we choose a two-level atom inside a photonic 
band gap (PBG) [21 BP] . Fictitious and pseudomode 
methods [T5J [T3] do not work for this system since the 
typical reservoir distribution function for PBG is not a 
meromorphic function due to the band edge. Moreover, 
an earlier attempt to develop a specific jump approach 
for this system [21] has been shown to be correct only 
in the Born-Markov limit [22j. One of the reasons for 
this is that the method of Ref. [5TJ fails to describe the 
reabsorption of photons by the atoms [55] . Our method 
succeeds in this by using non-Markovian quantum jumps 
[c.f. Eq. ([EJ]. This example also shows that local-in-time 
master equations can be used to describe non-Markovian 
dynamics for strong system-reservoir interactions. 

The master equation for the density matrix of the two- 
level system takes the form [T] 



m = 



ih 2 



A(t)(a-p(t)a + --{a + a-,p(t)} 



(8) 



where S(t) is the Lamb shift, A(t) the decay rate, er_ = 
|<?)(e|, and cr + = a_. Here, \g) denotes the ground state 



3 




5 10 
TIME [1/p] 



b £ 0.8 

" □ 

lu m 
< 
m 
O 

DC 
D_ 

LU 



< 

Q 



0.6 



0.4 



O 



0.2 



■ excited state 
ground state 



4 6 
TIME [1/p] 



10 



FIG. 1: (Color online) (a) The decay 
atom in photonic band gap as a function 
and exact results. In (a) the decay rate 
behavior with temporary negative values 
excited state probability of the atom and 
match between the exact and simulation 
pure states in examples 1 and 2 are \e) 
respectively. 



rate for a two-level 
of time, (b) NMQJ 
displays oscillatory 
In (b) we plot the 
the results show the 
results. The initial 
and {\g) + |e»/V2, 



FIG. 2: (Color online) An example realization with a jump 
- reverse jump cycle. The ground and excited state proba- 
bilities are given as function of time. The first jump at time 
t ~ 0.8 occurs at the positive decay rate region and destroys 
the superposition state. The second jump at t ~ 5.0 occurs 
at the negative decay region and recreates the superposition. 
The dotted lines show the evolution without any jumps. 



of the two-level atom, |e) the excited state, and there is 
one decay channel taking the atom from |e) to \g). We 
calculate the Lamb shift and the decay rates by using 
Eq. (2.21) of Ref. |D] and Eqs. (10.22) and (10.23) from 
Ref. [I]. The oscillatory behavior and negative values of 
the decay rate are displayed in Fig. [T] (a). 

Figure [T] (b) shows the match between the exact result 
[c.f. Eq. (2.21) of Ref. [20]] and the simulation results 
with N = 10 5 realizations for two different initial states. 
We have chosen parameters which correspond to Fig. 1 
of Ref. [5U] with the detuning 5 = —(3 from the edge of 
the gap. Here, (3 = (w^/ 2 d 2 /Qireohc 3 ) 2 ^ , where too is the 
Bohr frequency and d the absolute value of dipole mo- 
ment of the atom. The results illustrate a typical feature 
of PBG: atom-photon bound state and population trap- 
ping. Figure [2] displays an example of non-Markovian 
quantum jump in a single realization of the process for 
the case of initial superposition state. First, during the 
positive decay, a jump takes the atom to its ground state 
and the excitation resides in the environment. Later on 
with negative rate, the superposition state is restored by 
a non-Markovian quantum jump, and the photonic com- 
ponent is reabsorbed by the atom. 

Insight by NMQJ. In the PBG example above, the key 
ingredient to describe non-Markovian memory is the vir- 
tual photon emission-reabsorption cycle on the level of 
single realization. The physical state of the system is 
given by the density matrix, i.e., the ensemble of state 
vectors. This illustrates an interesting aspect of our 
method: it is possible to describe the effects of non- 
Markovian memory without extending the Hilbert space 
of the reduced system, which is a trait used in the pre- 
viously developed jump methods |T51 EE H3 [TH]. In 



NMQJ method, the memory of the ensemble member 
\ip a ), i-e., the information about the state before the posi- 
tive rate jump to the state \ip a ) occurred, is carried by the 
other ensemble member |Vv)- Consequently, the density 
matrix and the corresponding ensemble indeed carry in- 
formation on the earlier state of the system. 

Negative decay rates, which occur in non-Markovian 
systems, can be interpreted in the following way. Dur- 
ing the initial period of positive decay, the correspond- 
ing jumps distribute the state vector probability over the 
Hilbert space accordingly; the number of terms in the 
summation of Eq. ^ increases. When the decay rate 
later on becomes negative, which indicates the memory 
effects, the direction of the probability flow is reversed. 
This means that a process \ip) — > with negative rate 
corresponds to <— From classical perspective, 

it seems rather usual that changing the sign of the rate 
of the process means that the process occurs to the op- 
posite direction. In the quantum world with superpo- 
sitions, probability amplitudes and coherences the issue 
is less straightforward. In our method, this appears as 
a restoration of seemingly lost superpositions and subse- 
quent revival of coherences. 

The algorithm and numerical efficiency. Since in the 
NMQJ method the realizations depend on each other due 
to memory effects [c.f. Eq. pj], it seems at first sight that 
all the N ensemble members have to be evolved simul- 
taneously. However, according to Eq. the ensemble 
consists of several copies of each \tp a (t)}. Obviously, there 
is no need to have on a computer several copies of the 
same state vector. It is sufficient to have one copy and 
the corresponding integer number N a . Any number N 
of the realizations of the process can be done by mak- 
ing N c ff <§; N state vector evolutions where N e g is equal 
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to the number of terms in the summation N — J2 a ^a- 
When the realizations of the process are generated on 
a computer, a jump means changing the integer num- 
bers N a (t) accordingly in Eq. A considerable saving 
in CPU time is achieved since it is not necessary at each 
point of time evolve N state vectors, instead, it is enough 
to decide N times if the jumps occurred or not. 

Let us illustrate this with an example. In the PBG case 
above, we have N c g = 2 and the corresponding state vec- 
tors are: \ipo(t)) and \tpi(t)) = \g) for all t. These are the 
initial state affected by the deterministic evolution and 
the ground state, respectively. In the positive decay re- 
gion the jumps occur as |^o(*)) — * IV'iC*)) = eacn 
jump reduces No by 1 and increases N\ by 1. In the neg- 
ative decay region the process goes to opposite direction 
\ipi(t)) = \g) —> \tpo(t)}. In the optimized simulation to 
have 10 5 realizations we need to generate only one deter- 
ministic evolution for |^>o(*))> and then decide the jumps 
as described above. 

In QSD [12], the stochastic change of the state vectors 
is continuous which leads in practice to iV e flf = N. For the 
doubled Hilbert space (DHS) method [17], the norm of 
the state vectors increase in the negative decay region. As 
a consequence, the norm of a given state vector depends 
on the point of time where the DHS jump happens during 
the negative decay. In the ensemble, the jumps occur at 
each time point and N e g becomes large compared to the 
NMQJ method. Moreover, the DHS state vectors are 
evolved in the Hilbert space twice as large as in NMQJ. 

In contrast to the DHS method, the triple Hilbert space 
method (THS) preserves the norm of state vectors [15] . 
However, in the most general case, when the jump opera- 
tors depend on time in the master equation ([I]) , the jumps 
with the extended THS operators increase N e s at each 
point of time during the negative decay. Consequently, 
the THS method can not use the built-in optimization of 
the NMQJ method. Moreover, the THS method has two 
other ingredients which have an impact on its numeri- 
cal performance: (i) the need for 4 times larger number 
of decay channels than NMQJ uses [see the text below 
Eq. (57) in Ref. [15]]. (ii) the state vectors live in the 
space which is 3 times larger than the original one [see 
Eqs. (27)-(29) in Ref. [18J]. The consequent complica- 
tions of the THS method make it difficult to make a gen- 
eral statement on its numerical performance. However, 
all the facts above lead to the conclusion that even the 
most cautious estimate would give roughly an order of 
magnitude difference in the numerical efficiency between 
the NMQJ and the THS methods. 

Conclusions. The quantum jump description for 
Markovian systems (MCWF) is widely accepted due to 
its straightforward nature and the simple physical picture 
that it provides. For non-Markovian systems, the NMQJ 
method maps memory into reverse jumps that restore 
quantum superpositions. Furthermore, our approach be- 
comes equivalent to the standard MCWF method in 



the Markovian limit. In a broader view, the continu- 
ously growing interest in quantum information [6] and 
nanophysics [23] emphasises the need to consider single 
quantum systems at diminishing time scales and in tai- 
lored and finite environments. This development pro- 
vides the background for the NMQJ approach. 
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